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PREFACE 


I  have  always  been  interested  in  fluid  flows  and  ways 
of  analyzing  them.  When  I  was  offered  the  opportunity  by 
Maj  Hodge  to  investigate  a  method  of  analyzing  the 
noni sothermal  wall  effect  in  hypersonic  flow,  I  was  both 
excited  and  interested  in  the  project.  The  study  itself  was 
intended  to  be  a  proof-of-concept  for  the  numerical  methods 
used  and  a  stage  to  show  the  aerospace  industry  that  the 
nonisothermal  wall  effect  is  important  and  can  be  included 
in  design  analyses.  While  the  results  of  the  thesis  have 
shown  that  the  particular  method  used  is  not  the  optimal 
approach,  I  feel  that  the  objectives  of  the  study  have  been 
met . 

Attempting  any  project  of  this  size  creates  a  long  list 
of  indebtedness;  the  fact  that  is  was  necessary  for  me  to 
continue  the  effort  after  leaving  AFIT  only  lengthens  that 
list.  One  of  the  most  important  people,  of  course,  has  been 
my  advisor,  Maj  James  K.  Hodge.  His  insight,  direction,  and 
especially  patience  have  been  invaluable  to  me  in  completing 
this  project.  The  help  I  received  from  Dr  (Jrmilla  Ghia  of 
the  University  of  Cincinnati  while  she  was  at  AFIT  as  a 
Distinguished  Visiting  Professor,  and  from  Dr  Wilbur  Hankey 
and  Dr  Joseph  Shang  of  the  Flight  Dynamics  Laboratory  at 
Wright-Patterson  AFB  has  been  valuable  and  illuminating.  I 
thank  Dr  Ghia  for  introducing  me  to  the  rigor  necessary  to 


wrest  correct  results  from  recalcitrant  equations.  Drs 
Hankey  and  Shang  were  invaluable  in  helping  master  the 
nuances  of  computational  fluid  dynamics. 

Finishing  a  thesis  after  leaving  AFIT  requires  the 
active  support  from  the  people  you  work  with.  First  on  this 
list  is  my  branch  chief,  Mr  Ed  Barth.  Without  his  active 
support  of  my  studies,  this  thesis  would  not  have  been 
completed.  Drs  Phil  Kessel  and  Joseph  Baum  have  been  very 
helpful  in  aiding  and  critiquing  my  efforts.  Finally,  the 
people  in  my  section,  who  understood  when  their  section 
chief  was  somewhat  preoccupied  at  times. 

The  most  important  acknowledgment  has  been  saved  for 
last.  It  is  an  oft-stated  truism  that  "I  couldn't  have  done 
it  without  my  wife."  This  is  a  truth  that  repetition  cannot 
dim.  Mary  has  been  at  my  side  throughout  AFIT  and  after, 
egging  me  on,  supporting  me,  and  showing  me  that  someone 
else  cared,  too.  More  than  anyone  else,  Mary,  this  one’s 
for  you.  And  to  David,  who  shouldn't  have  but  did. 


Timothy  K.  Roberts 
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ABSTRACT 


Wind  tunnel  tests  of  Space  Shuttle  Orbiter  insulating 
articles  have  demonstrated  the  presence  of  a  nonisothermal 
wall  effect,  which  is  a  lag  in  heat  transfer  recovery  after 
the  flow  passes  over  a  surface  temperature  discontinuity 
resulting  in  a  downstream  transport  of  energy.  Theoretical 
analyses  and  numerical  simulations  of  hypersonic  flow  over 
discontinuous  nonisothermal  surfaces  using  boundary  layer 
theory  have  also  indicated  the  presence  of  this  effect. 

This  thesis  studies  the  nonisothermal  wall  effect  by 
modeling  the  hypersonic  flow  over  an  inclined  wedge  with  a 
discontinuous  nonisothermal  surface.  The  flow  is  modeled 
using  the  two-dimensional  Navier-Stok.es  equations.  MacCor- 
mack's  method  is  used  to  solve  the  Navier-Stokes  equations. 
The. program  used  to  implement  these  methodologies  is  discus¬ 
sed  and  a  listing  given.  A  sem  i-adapt  i  ve  grid  is  used  to 
represent  the  physical  conditions  of  the  problem.  Heat 
transfer  is  presented  as  a  nondimensional  ratio  of  the  local 
convective  heat  transfer  coefficient  to  a  reference  heat 
transfer  coefficient.  <  •  ,  -  ,  <  < 

The  results  of  this  study  show  that  the  nonisothermal 
wall  effect  cn  be  successfully  modeled  using  the  two-dimen¬ 
sional  Navier-Stokes  equations  and  MacCormack's  explicit 
method.  The  lag  in  the  recovery  of  the  convective  heat 
transfer  coefficient  is  found  to  match  lags  seen  in  other 


analyses  of  the  problem.  Due  to  the  high  Mach  number 
modeled,  a  shock  wave-boundary  layer  interaction  is  found  to 
have  an  effect  on  the  heat  transfer.  A  significant  surface 
pressure  spike  is  found  to  occur  downstream  of  the 
discontinuity,  which  may  be  a  result  of  the  increase  in  size 
of  the  momentum  and  thermal  boundary  layers  at  the 
discontinuity. 

The  study  concludes  that  the  nonisothermal  wall  effect 
can  be  adequately  modeled  by  the  two-dimensional  Navier- 
Stokes  equations;  that  the  shock  wave-boundary  layer 
interaction  does  have  an  effect  on  the  heat  transfer;  and 
that  the  occurence  of  a  spike  in  surface  pressure  may  be  a 
unique  result  of  the  nonisothermal  wall  effect.  Significant 
resommendat  ions  include  the  need  for  further  study  of  the 
nonisothermal  wall  effect,  the  need  to  use  more  optimal 
grids  and  solution  methods,  the  need  to  more  thoroughly 
investigate  the  shock  wave-boundary  layer  effect,  and  the 
need  for  further  study  of  the  surface  pressure  response  to 


the  nonisothermal  wall  effect. 


A  NUMERICAL  SOLUTION  OP  A  NONISOTHERMAL  WALL 


USING  THE  TWO-DIMENSIONAL  NAVI ER-STOKES  EQUATIONS 


CHAPTER  I :  INTRODUCTION 


BACKGROUND 

The  development  and  use  of  the  Space  Transportation 
System  has  given  rise  to  some  unanticipated  problems  in 
reentry  aerodynamics.  The  fact  that  the  flow  field  over  the 
Orbiter  is  laminar  during  reentry  down  to  Mach  8  permits 
some  phenomena  usually  masked  by  turbulent  flow  to  become 
evident.  One  such  phenomenon  which  significantly  affects 
the  flow  and  the  heat  transfer  is  the  nonisothermal  wall 
effect.  This  effect  occurs  when  flow  passes  over  a  surface 
composed  of  two  or  more  materials  of  different  thermal 
response.  The  nonisothermal  wall  effect  occurs  at  numerous 
locations  on  the  Orbiter's  surface.  One  such  location  is 
the  junction  between  the  reusable  carbon-carbon  (RCC)  nose- 
cap  and  the  reusable  surface  insulation  (RSI)  tiles.  The 
flow  must  pass  over  a  discontinuity  in  wall  temperature  due 
to  different  thermal  responses.  Similar  discontinuities 
exist  at  numerous  other  places  on  the  Orbiter's  surface. 

During  reentries  of  earlier  space  vehicles  such  as 
Apollo  or  ballistic  missile  reentry  vehicles,  the  flow  over 
the  vehicle  was  turbulent.  The  mixing  inherent  in  turbulent 
flow,  augmented  by  the  surface  roughness  of  these  vehicles, 

) 


reduces  the  dramatic  change  in  the  wall  heat  transfer  in¬ 
duced  by  a  nonisothermal  wall  condition.  However,  an  Orbi- 
ter  reentry  is  characterized  by  a  laminar  flow  field  over 
much  of  the  vehicle's  surface  until  the  Orbiter  has 
decelerated  to  about  Mach  8;  at  that  point,  transition  to 
turbulent  flow  is  abrupt.  Thus,  the  nonisothermal  wall 
effect  becomes  evident.  In  fact,  the  localized  scorching 
and  discoloration  seen  on  Columbia,  Challenger,  and  Discove¬ 
ry  seem  to  be  caused  by  the  nonisothermal  wall  effect. 

Studies  of  wind  tunnel  tests  of  Orbiter  insulating 
articles  have  demonstrated  the  presence  of  this  effect. 
Theoretical  analyses  using  boundary  layer  theory  of 
hypersonic  nonisothermal  walls  had  indicated  that  the 
discontinuity  in  wall  temperature  would  induce  changes  in 
the  flow.  Numerical  simulations  using  classical  boundary- 
layer  theory  backed  up  these  findings.  However,  the  wind 
tunnel  tests  indicated  that  when  the  flow  passed  from  a 
steel  surface  to  a  surface  of  Orbiter  insulation  material, 
wall  temperature  recovery  was  much  slower  than  that 
predicted  by  either  Eckert  high-speed  plate  theory  or  boun¬ 
dary-layer  simulations.  The  data  gathered  in  these  tests 
would  seem  to  indicate  that  some  other  effect  is  at  least  as 
important  as  normal  diffusion  in  transporting  energy. 
First,  axial  diffusion  of  energy  at  a  location  where  there 
is  a  discrete  change  in  wall  temperature  over  an  axial 
distance  of  order  less  than  the  normal  boundary  layer 


thickness  is  one  possible  effect.  Second,  the  viscous  shock 
wave  -  boundary  layer  interaction  at  hypersonic  speeds 
changes  the  pressure  and  heat  transfer  in  an  axial  direc¬ 
tion.  These  two  effects  have  not  been  investigated  on  a 
nonisothermal  wall. 

OBJECTIVE 

The  objective  of  this  thesis  is  to  numerically  model 
hypersonic  flow  over  an  isothermal  and  a  nonisothermal  wall 
using  the  two-dimensional  Navier-Stokes  equations.  The  use 
of  the  Navier-Stokes  equations  allows  accurate  modeling  of 
the  axial  diffusion  if  there  is  sufficient  axial  grid  reso¬ 
lution.  The  nonisothermal  wall  effect  should  be  very 
noticeable.  This  paper  will  attempt  to  show  the  effects  of 
the  shock  wave-boundary  layer  interaction  and  axial  diffu¬ 
sion  of  energy  on  the  nonisothermal  wall  effect. 


CHAPTER  II:  THEORETICAL  DEVELOPMENT 


OVERVIEW 

The  objective  of  this  thesis  is  to  numerically  model 
hypersonic  flow  over  a  wedge  with  and  without  a  discontinui¬ 
ty  in  surface  temperature.  This  chapter  will  discuss  the 
theory  behind  the  investigation.  The  problem  will  be 
described  and  possible  avenues  of  investigation  proposed. 
The  best  method  of  solution  will  be  selected  and  reasons  for 
that  selection  will  be  discussed. 

PROBLEM  DESCRIPTION 

The  physical  problem  to  be  modeled  is  a  wedge  in  a 
hypersonic  flow  field.  The  problem  is  described  in  Cappela- 
no  (1:2-5),  Hodge  et  al  (2:2),  and  Woo  (3).  See  Figure  1 
for  a  picture.  A  steel  test  article  has  a  6  inch  by  6  inch 
test  plate  inserted  in  it  7  inches  behind  the  leading  edge 
of  the  test  article.  Cappelano  reviewed  and  interpreted 
wind  tunnel  test  data  for  three  types  of  inserts:  a  thin 
steel  plate  (for  calibration),  a  piece  of  high  temperature 
reusable  surface  insulation  (HRSI),  and  a  piece  of  flexible 
reusable  surface  insulation  (FRSI).  Woo  reduced  the  data 
for  the  HRSI  and  the  FRSI  test  runs.  Each  insert  was  in¬ 
strumented  with  varying  numbers  of  surface  thermocouples  to 
record  the  axial  surface  temperature  distribution.  The  FRSI 
insert  had  the  most  thermocouples  and  hence  the  most  com¬ 
plete  set  of  axial  surface  temperature  readings;  thus,  the 
FRSI  insert  is  modeled  in  this  study.  Only  one  angle  of 
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attack  is  simulated  in  this  study;  restrictions  on  computer 
time  required  that  excursions  to  other  angles  of  attack  be 
deferred  to  later  studies. 


The  flow  conditions  modeled  were  as  follows:  a 
freestream  Mach  number  of  14.24,  freestream  stagnation  tem¬ 
perature  of  2000  degrees  Rankine  (deg  R) ,  and  a  wedge  incli¬ 
nation  angle  of  3  degrees.  The  combination  of  high  Mach 
number  and  stagnation  temperature  produce  a  very  low 
freestream  static  temperature  of  48  deg  R;  several  problems 
will  arise  from  this.  One  of  the  immediate  ones  is  that 
gradients  will  be  very  steep  in  the  boundary  layer, 
requiring  significant  resolution  over  a  very  thin  layer. 
Another  problem  is  that  any  numerically-induced  waves  could 
easily  cause  temperatures  just  in  front  of  the  leading  edge 
to  become  negative  in  the  computations,  thus  causing  the 
solution  to  stop. 

AVENUES  OF  INVESTIGATION 

Several  avenues  of  investigation  are  available  to  solve 
a  problem  such  as  this,  depending  on  the  desired  results. 
Approaches  fall  into  two  broad  categories:  boundary  layer 
solutions  and  Navier-Stokes  solutions.  Both  approaches  will 
be  examined  in  some  detail  below.  Boundary-layer  oriented 
solutions  offer  a  relatively  fast  solution,  but  at  the 
expense  of  some  details  of  the  flow.  If  the  problem  is  such 
that  those  approximations  cannot  be  accepted,  then  the 
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investigator  must  use  methods  of  increasing  complexity, 
collectively  known  as  the  Navier-Stokes  family  of  solutions. 

There  are  several  levels  of  complexity  within  the  fami¬ 
ly  of  Navier-Stokes  solutions.  Slightly  more  complicated, 
but  consequently  more  realistic,  are  the  parabolized  Navier- 
Stokes  equations.  Beyond  that  are  the  "full"  Navier-Stokes 
equations  in  two  or  three  dimensions.  For  reasons  to  be 
explained  in  more  detail  below,  a  two-dimensional  Navier- 
Stokes  solution  was  selected  for  this  problem. 


BOUNDARY  LAYER  SOLUTIONS 

Boundary  layer  solutions  depend  on  Prandtl's  boundary 
layer  theory.  The  approximations  made  by  this  theory  are 
such  that,  in  most  cases,  the  solution  resulting  from  boun¬ 
dary  layer  theory  is  very  close  to  the  actual  flow  field. 
This  is  especially  true  in  low  Mach  number  or  incompressible 
flows.  At  hypersonic  Mach  numbers,  the  boundary  layer  ap¬ 
proximations  tend  to  give  results  that  diverge  somewhat  from 
the  actual  flow  field,  especially  when  there  is  an  interac¬ 
tion  between  the  oblique  shock  and  the  viscous  boundary 
layer . 

The  boundary  layer  equations  developed  by  Prandtl  make 
several  important  assumptions.  These  assumptions  follow 
from  some  basic  observations  about  the  nature  of  any  flow 
field.  The  boundary  layer  is  defined  as  that  region  near 
the  surface  of  the  body  in  which  the  flow  velocity  increases 


from  zero  at  the  wall  (due  to  the  viscous  nature  of  fluids) 


to  the  freestream  velocity  at  some  distance  from  the  body. 
Experimentation  has  shown  that  this  distance,  the  boundary 
layer  thickness,  is  normally  quite  small  compared  to  other 
character istic  lengths  of  the  problem.  Prandtl's  boun¬ 

dary  layer  equations  can  be  derived  from  the  full  set  of  the 
two-dimensional  Navier-Stokes  equations  using  order  of  mag¬ 
nitude  analysis.  The  variables  used  will  be  assumed  to  be 
in  nond imens i ona 1  form.  For  example,  the  axial  velocity 
component  is  nond i me  ns  i  ona 1 i ze d  with  respect  to  the 
freestream  velocity  and  has  an  order  of  magnitude  of  one. 
Likewise,  the  axial  length  is  nond imens ional i zed  with 
respect  to  the  length  of  the  wedge  and  is  of  order  one.  The 
normal  length  is  nond imens iona 1 ized  with  respect  to  the 
boundary  layer  thickness,  delta.  Thus,  it  has  an  order  of 
magnitude  of  delta.  With  these  building  blocks,  the  two- 
dimensional  Navier-Stokes  equations  can  be  reduced  to 
Prandtl's  boundary  layer  equations. 

Consider  first  the  continuity  equation: 

4  r  0  (2.1) 

l  « 

The  appropriate  order  of  magnitude  has  been  written  below 
each  term.  In  order  to  maintain  the  continuity  equation,  it 
can  be  seen  that  the  order  of  magnitude  of  the  normal  velo¬ 
city  component,  v,  must  be  delta. 

Next,  consider  the  axial  linear  momentum  equation: 
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2.2 


The  stress  tensor  terms  cPq  and  are  defined  as  (4:61): 
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Substituting  these  terms  into  Eq.  2.2: 
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The  appropriate  order  of  magnitude  has  been  written  under 
each  term  for  which  it  is  known.  By  discarding  those  terms 
with  orders  of  less  than  one  (i.e.,  orders  of  delta  and  its 

multiples),  the  axial  linear  momentum  equation  becomes: 
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Now  consider  the  normal  linear  momentum  equation: 
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The  new  stress  tensor  term  i-s  defined  as  (4:61): 
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Substituting  for  the  stress  tensor  terms  yields: 
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where  the  orders  of  magnitude  were  again  written  under  each 


term,  as  usual.  The  only  term  to  survive  in  this  equation 
is  the  partial  derivative  of  pressure  in  the  y-direction: 
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Finally,  consider  the  energy  equation: 
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Substituting  once  more  for  the  stress  tensor  terms: 
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i  i  a  i  '  a  ,  c 


(2.10) 


"  v  (ile  )i]  4  ^  [/°ve  ‘  Pf  U  ~  u  («i«C 

s  ^  i  \  i  ^  r  i  sl 


i 

/  ■)  tr  „  ^y 

e  c)> 

i  i  i 

&  < 


(2.11) 


the  order  of  magnitude  are  once  again  listed  under  each  term 
for  which  they  are  known.  The  resulting  equation  is: 

h  „  j£  ,  v»s  :  .  _J —  .  J=L 

Ot  Ox  ^  P  Z+  J  (2.12) 

This  is  equivalent  to  the  usually  accepted  boundary  layer 
energy  equation. 

Eqs.  2.1,  2.5,  2.9,  and  2.12  show  explicitly  the  major 
assumptions  in  Prandtl's  boundary  layer  theory.  Eq.  2.1 
shows  that  the  dominant  velocity  term  is  in  the  axial  direc¬ 
tion;  it  will  dominate  the  characteristics  of  the  boundary 
layer.  Eq.  2.5  reinforces  this  point  and,  along  with  Eq. 
2.9,  shows  that  the  pressure  distribution  in  the  boundary 
layer  varies  only  in  an  axial  direction.  The  surface  pres¬ 
sure  is  virtually  equal  to  the  impressed  freestream  pres¬ 


sure.  Any  shock  wave-boundary  layer  interaction  is  ignored 
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by  this  assumption.  Finally,  Eq.  2.12  shows  that  virtually 
all  the  heat  transfer  occurs  in  the  normal  direction,  with 
almost  no  axial  diffusion.  The  consequence  5  of  Eq.  2.12  are 
that  ax  is  of  order  one.  But  if  Ax  is  of  order  delta,  such 
as  across  a  discontinuity  in  wall  temperature,  the  axial 
diffusion  of  energy  is  neglected.  It  should  be  noted  that 
these  results  are  best  applied  to  incompressible  or  low- 
speed  compressible  flow  problems. 

REASONS  FOR  DISCARDING  BOUNDARY  LAYER  SOLUTIONS 

Boundary-layer  oriented  solutions  were  discarded  a 
priori  for  this  study  for  two  basic  reasons.  First,  the 
temperature  derivatives  in  the  streamwise  or  axial  direction 
are  ignored.  Second,  the  boundary  layer  solutions  do  not 
adequately  model  the  shock  wave-boundary  layer  interaction 
that  occurs  at  hypersonic  Mach  numbers. 

The  question  of  temperature  derivatives  is  straightfor¬ 
ward.  Recall  that  this  study  wishes  to  investigate  the 
nonisothermal  wall  effect.  One  of  the  constituents  of  the 
non i sothe rma 1  wall  effect  could  be  axial  diffusion.  Yet 
boundary  layer  theory  ignores  axial  diffusion.  Therefore, 
nonisothermal  wall  effect  cannot  be  properly  investigated 
with  a  tool  that  ignores  axial  diffusion. 

The  interaction  of  the  leading  edge  shock  wave  and  the 
boundary  layer  is  a  important  part  of  the  the  flow  field 
structure.  Such  shock  wave-boundary  layer  interactions  have 
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long  been  recognized  as  important;  see,  for  example,  an 
early  article  by  Nagamatsu  and  Sheer  (  4  :4  54-462).  In  this 
article,  Nagamatsu  described  how  induced  surface  pressures 
from  the  interaction  were  as  high  as  twelve  times  the 
freestream  pressure  for  a  similar  freestream  Mach  number. 
Clearly,  in  such  an  environment,  boundary  layer  assumptions 
break  down  and  more  realistic  models  must  be  used.  These 
two  factors  were  the  greatest  contributors  to  discarding 
boundary  layer  solutions  for  this  effort. 


NAVI ER-STOKES  SOLUTIONS 

The  method  of  choice  for  solving  fluid  flow  problems  is 
solving  the  Nav ie r-S tokes  equations.  Unfortunately,  a 
general  analytical  solution  is  not  available  for  the 
majority  of  applications;  numerical  solutions  must  be 
resorted  to,  as  they  were  in  this  study.  'Small-scale' 
effects  such  as  the  non  i  so  the  rma  1  wall  effect  can  be  seen 
and  their  importance  evaluated. 

The  Na v ie r-S toke s  equations  used  here  are  written  in 
conservative  form,  that  is,  in  terms  of  p ,  p^<-  ,  p\i  , 
and  /3d.  The  full  set  of  equations  used  are: 


r  - 


/ouv  -  T(. 


I  <}T 


TL  ~  *  vr„  -  v  r . 


where : 


-  u  -  vr 


=  -  ?  *  1  /  ^  ^ 


(2.13a) 


_  _  _  P  _L  /  H  ^  _  £  iv  \ 

°u  1  *  Jit  k  3  d*  3  ,)x,  J 

.  v  x  (2.13a) 

V  .  J_  /£*  *  h  ) 

^  Rri.  v  <?7  <3*  / 

-V.  *  /  M.  <)v  _  2.  e)  u.  \ 

'  ~r*  (ieUav,  1z>  ) 

The  grid  used  in  this  study  is  semi-adaptive;  a  complete 
discussion  will  be  given  below.  The  code  used  tranforms  the 
equations  internally. 


HEAT  TRANSFER  MODELING 

Heat  transfer  is  expressed  in  this  study  as  the  convec¬ 
tive  heat  transfer  coefficient,  h,  also  known  as  the  film 


and  is  expressed  in  B t u/ f t ^ - se c-deg  R.  To  be  able  to  com¬ 
pare  results  with  Woo  (3)  and  Hodge,  the  convective  heat 
transfer  coefficient  is  nondimensionalized  with  respect  to  a 


reference  convective  heat  transfer  coefficient  based  on 


free stream  quantities: 


kft(  --  C.13Z  Op 


^  poo  '-fjo  , 

*  / 


‘/z. 


(2.15) 


Note  that  the  only  variable  is  the  distance  along  the  wedge 
in  feet,  x.  Thus,  hre^  is  truly  a  reference  condition  at 
any  point  along  the  surface  of  the  wedge.  According  to 
Eckert,  for  an  isothermal  wall  and  no  viscous  shock-boundary 
layer  interaction,  the  ratio  of  h  to  href  is  independent  of 
x . 


In  both  quantities,  the  coefficient  of  viscosity,  is 

determined  by  Sutherland's  law: 

_  Vi. 


2 .2-7  >.(  o 


(2.16) 


Although  Sutherland's  law  is  generally  considered  to  be 
valid  down  to  about  180  deg  R  (4:19),  Hodge  et  al  have 
determined  that  Sutherland's  law  can  be  applied  to  lower 
temperature  flows  with  minimal  error  (2:2).  The  thermal 


conductivity  of  air  is  computed  by: 


J  -  -in 


(2.17) 


JT  ?r  '  ti+u. 

with  units  of  Btu/f t-sec-deg  R.  The  Prandtl  number  is  taken 

to  be  constant,  Pr  =  0.72.  Both  the  coefficient  of  viscosi¬ 
ty  and  thermal  conductivity  are  computed  at  the  wall  condi¬ 
tions  . 

* 

Eckert's  reference  temperature,  T  ,  can  be  used  to 
validate  the  use  of  Sutherland's  law  for  computation  of  the 


coefficent  of  viscosity.  T  is  defined  as  (6:304): 

T'-  (2.18) 

where  Te  is  the  boundary  layer  edge  temperature  (also  equal 
to  the  freestream  static  temperature  behind  the  shock),  and 
T  w  is  the  wall  temperature  (which  will  vary  for  the 
non i so t he rm a  1  wall).  Adiabatic  wall  temperature,  Taw,  is 
defined  as  (6:301): 

i  (2.19) 

where  the  recovery  factor,  r,  is  taken  to  be  a  constant,  r  = 
0.9.  This  is  not  the  generally  accepted  definition;  the 
recovery  factor  is  usually  defined  as  the  square  root  of  the 
Prandtl  number.  However,  this  alternative  definition  sim¬ 
plifies  computations,  causes  only  a  few  percent  deviation 
from  convective  heat  transfer  coefficient  computed  with  the 

"traditional"  recovery  factor,  and  allows  direct  comparison 

* 

with  Hodge's  and  Woo's  results.  T  is  an  empirical 

correcting  factor.  Eckert  found  that  if  one  assumes  the 

* 

specific  heats  to  be  constant,  one  can  use  T  to  correlate 

many  "exact"  laminar  boundary  layer  solutions  to  within  a 

* 

few  percent.  It  is  also  apparent  that  if  T  is  greater  than 

180  deg  R,  Sutherland's  law  is  appropriate  for  the  coputa- 

* 

tion  of  the  coefficient  of  viscosity.  In  this  case,  T  is 
equal  to  over  765  deg  R.  Thus,  the  use  of  Sutherland's  law 
is  justified. 
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CHAPTER  III:  NUMERICAL  SOLUTION  METHODS 


OVERVIEW 

This  chapter  will  consider  the  generation  and  use  of 
the  numerical  grid  and  the  numerical  implementation  of  the 
Navier-Stokes  equations  in  NONISOCODE  using  MacCormack's 
explicit  method. 


NUMERICAL  GRID 

The  numerical  grid  used  in  this  study  was  a  semi- 
adaptive  grid.  It  was  generated  in  such  a  way  so  as  to 
achieve  two  goals:  a)  maintain  about  ten  points  in  the 

boundary  layer  at  all  positions  along  the  plate,  and  b) 
maintain  thirty  points  in  the  normal  direction  at  all  sta¬ 
tions  along  the  plate.  The  combination  of  these  two  goals 
led  to  the  use  of  a  semi-adaptive  grid. 

The  normal  grid  distribution  was  found  to  be  the  solu¬ 


tion  of  the 


The  solution 


following  differential 

&  - -  ° 
to  this  equation  is: 


equation: 


(3.1) 


(3.2) 


where  y  is  the  solved  physical  location  of  the  normal  grid 
point,  ymax  is  the  highest  physical  grid  location,  ^  is 
athe  number  of  the  grid  point  being  solved  for,  max  is 
the  maximum  number  of  grid  points,  and  Q  is  the  stretching 


exponent . 

The  primary  goal  was  to  maintain  about  ten  points  in 
the  boundary  layer  at  all  stations  in  the  streamwise  direc¬ 
tion.  This  was  required  in  order  to  adequately  model  the 
temperature  inversion  that  exists  there.  An  earlier  phase 
of  this  effort  had  shown  that  it  is  very  easy  to  miss  that 
inversion  if  the  grid  isn't  carefully  planned,  thus 
rendering  the  results  useless.  The  limit  of  ten  points  in 
the  boundary  layer  was  set  to  keep  the  required  stretching 
down;  NONISOCODE  will  not  run  a  grid  that  has  a  stretching 
exponent  much  greater  than  -0.15.  Ten  points  was  considered 
an  acceptable  compromise  between  desired  level  of  detail  and 
the  requirements  of  the  parent  code. 

The  secondary  goal  of  maintaining  thirty  points  in  the 
normal  direction  at  all  stations  was  an  objective  of 
convenience.  It  would  have  been  possible  to  use  forty, 
fifty,  or  more  points  in  the  normal  direction  to  achieve  any 
desired  level  of  detail;  however,  the  price  paid  is  that  of 
increased  run  time  for  the  solution  and  increasingly 
unwieldy  arrays.  Use  of  more  than  thirty  points  in  the 
normal  direction  was  considered  unneccesary  for  the  level  of 
detail  desired.  To  do  so  would  have  required  more  modifica¬ 
tions  to  the  code;  simple  enough  in  theory  but  prone  to 
error  in  application. 

Generating  the  grid  was  a  trial-and-error  process. 
There  were  several  known  constraints.  One  was  the  minimum 


height  of  the  grid  needed  at  the  downstream  edge  of  the 

computational  domain  in  order  to  pass  shock  and  Mach  waves 

(i.e.,  not  have  them  reflect  off  the  far-field  boundary). 

Another  constraint  was  the  compromise  between  the  number  of 

points  in  the  boundary  layer  and  the  total  number  of  points 

in  the  normal  direction  mentioned  above.  It  was  decided 

that  the  grid  should  grow  in  the  normal  direction 

proportional  to  the  boundary  layer  thickness.  Allowing  ymax 

to  increase  while  maintaining  the  same  exponential  stretch 

and  total  number  of  points  was  done  using  this  equation: 

.  Vj- 


GRIDY i  =  GRIDYref 


(3.3) 


The  normal  coordinate  is  proportional  to  the  square  root  of 
the  ratio  of  the  current  streamwise  coordinate  to  a 
reference  coordinate.  The  constant  of  proportionality  is 
the  reference  normal  coordinate.  The  choice  of  the 
reference  coordinate  is  somewhat  arbitrary;  the  driving 
factor  influencing  it  is  the  location  of  the  closest  point 
to  the  leading  edge  that  allows  ten  points  in  the  tempera¬ 
ture  inversion.  For  an  exponential  stretching  coefficient 
of  0  =  -0.15,  that  location  was  X/L  =  0.  The  actual  grid 
growth  began  at  the  leading  edge;  prior  to  that,  the  grid 
used  was  a  rectangular,  regular  grid  for  the  three  points  in 
the  freestream  ahead  of  the  wedge.  The  rationale  here  is 
that  freestream  conditions  are  not  computationally  severe 


and  hence  can  be  accomodated  in  a  regular  grid.  Beginning 


the  semi-adapt ive  grid  at  the  leading  edge  accounts  for  the 
growth  of  the  boundary  layer. 


Figure  2  is  a  plot  of  the  grid  used  in  this  study. 
Note  that  the  initial  gradient  of  the  parabolic  grid  is 
quite  steep.  This  is  because  the  boundary  layer  is  growing 
very  quickly  near  the  leading  edge.  Note  also  that  the 
steepest  grid  gradients  are  in  the  most  benign  flow  region. 
Near  the  leading  edge,  the  shock  wave  is  very  close  to  the 
surface  and  the  boundary  layer  is  very  thin.  Thus,  at  the 
far-field  boundary,  the  flow  is  essentially  inviscid  and  at 
freestream  conditions.  The  stagnation  region,  the  area  of 
sharp  turning  in  the  flow,  and  the  shock  wave  are  all  pre¬ 
sent  in  a  region  of  relatively  small  grid  gradients,  not 
unlike  the  gradients  found  elsewhere  in  the  computational 
domain . 


MACCORMACK'S  EXPLICIT  METHOD 

NONISOCODE  implements  the  Navier-Stokes  equations  using 
MacCormack's  explicit  method.  MacCormack's  method  is  one  of 
the  best  known  methods  for  numerically  solving  supersonic 
and  hypersonic  flow  problems.  Originally  presented  in  1969 
(8),  the  method  has  become  widely  used  for  many  high-speed 
problems . 

The  adjective  ’explicit*  refers  to  the  fact  that  a 
solution  at  any  point  in  the  flow  is  independent  of  the 

solution  at  any  other  point  at  that  time  step.  One  of  the 


distinguishing  characteristics  of  super-  and  hypersonic 
problems  is  that  information  cannot  travel  upstream,  since 
the  flow  velocity  everywhere  (except  in  a  very  thin  region 
within  the  boundary  layer)  is  greater  than  the  speed  at 
which  flow  information  can  travel  -  the  speed  of  sound.  in 
the  problem  under  consideration,  the  vast  majority  of  the 
flow  to  be  solved  is  in  the  hypersonic  flow  regime  with  only 
a  small  (but  important)  region  in  the  subsonic  regime. 
Thus,  a  solution  method  that  is  efficient  in  the  hypersonic 
regime  is  desirable.  MacCormack's  explicit  method  is  an 
efficient  method  at  high  velocities  (i.e.,  super-  and  hyper¬ 
sonic).  It  is  true  that  a  very  important  portion  of  the 
flow  to  be  solved  is  subsonic  (the  boundary  layer  where  the 
bulk  of  the  heat  transfer  occurs).  It  is  also  true  that 
implicit  methods  are  more  efficient  in  subsonic  regimes. 
However,  since  the  overall  flow  is  mixed,  MacCormack's  ex¬ 
plicit  method  is  the  more  advantageous  overall. 

DERIVATION  OF  MACCORMACK'S  EXPLICIT  METHOD 

Recall  the  two-dimensional  Nav ier-Stokes  equations 


derived  earlier: 


where : 
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(3.4a) 
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MacCormack's  method  states  that,  if  the  solution  Un<  is 

1  •  J 


(3.4c) 


known  at  time  t=n  t  at  each  point  in  the  mesh,  then  at 
time  t=(n+l)A,t  the  solution  can  be  described  as: 


U 


n+1,  .  =  *  (  £>  t)0n. 
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where  ^  (  At)  defines  these  operations: 


(3.5) 
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An  example  of  the  application  of  MacCormack's  explicit 


method  is  shown  by  considering  the  computation  of  density: 

/C-  *-  (3’7’ 

Now  we  write  Equation  3.7  in  terms  of  its  components: 
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Note  that: 
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(3.9) 

Thus,  the  two  computations  should  proceed  simultaneously. 


c* 


TRUNCATION  ERROR  OF  MACCORMACK ' S  EXPLICIT  METHOD 

The  truncation  error  of  MacCormack's  explicit  method  is 
second  order,  as  mentioned  earlier.  This  is  true  only  of 
the  combined  form: 

.  n*i  I  T  n  .  ,  r\  At  ,c  n  \ 

Uc-  -  ^  j  )  ‘  ^  j  “  »)  (3  10) 

J  Jt/H  J  V  ir.ri.  J  u-10) 

"  F-j  j  “  Sy  tKj,.  *V.  )\ 

By  themselves ,  the  predictor  and  corrector  are  only  first 
order  accurate.  The  use  of  backward  space  differencing  in 
the  predictor  and  forward  space  differencing  in  the  correc¬ 
tor  (or  vice  versa;  as  long  as  the  directions  are  different) 
allow  the  method  to  achieve  second  order  accuracy. 


STABILITY  OF  MACCORMACK'S  EXPLICIT  METHOD 

The  stability  of  MacCormack's  explicit  method  has  not 
been  analytically  demonstrated  yet.  According  to  MacCormack 
((9:3)  and  (9:152-154)),  some  insight  could  be  derived  from 
approximations.  For  instance,  one  could  linearize  the  Na- 
vier-Stokes  equations  and  then  study  the  amplification  of 
Fourier  components  of  the  solution  (9:152).  Alternatively, 
one  could  examine  three  separate  parts  of  the  Navier-Stokes 

equations  (inviscid,  viscous,  and  mixed  derivatives)  and 
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find  stability  conditions  for  each  part  (9:3).  MacCormack 
did  postulate  stability  criteria  as  follows: 


6  A 


&t  - 
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) 


( 3.11  ) 


v  Pr  &y 

This  condition  is  usually  called  the  Courant-Fr iedr ich-Levy 
( CFL)  condition.  The  implementation  of  MacCormack's  method 
in  NONISOCODE  uses  a  form  of  the  CFL  condition  computed  at 
each  streamwise  strip  for  different  Ay's  which  ignores  the 
viscous  contributions.  The  implementation  in  NONISOCODE  is 
also  in  transformed  coordinates,  so  that  its  form  is 
somewhat  different  from  that  of  Eq.  3.11. 
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CHAPTER  IV:  PROGRAM  DESCRIPTION 

OVERVIEW 

This  chapter  will  discuss  the  t'-’~  primary  codes  used  in 
this  study,  NONISOCODE  and  STEADY.  Each  code  will  be 
discussed  in  some  detail.  Listings  of  both  programs  can  be 
found  in  the  Appendix  A. 

NONISOCODE 

The  Nav  ier-Stokes  simulation  code  used  in  this  study, 
NONISOCODE,  is  a  derivative  of  a  code  written  by  Dr  J.  S. 
Shang  et  al  at  the  Flight  Dynamics  Laboratory  ( AFWAL/FI M ) , 
Wright-Patterson  AFB ,  Ohio.  A  copy  of  the  program  listing 
is  found  in  Appendix  A.  This  section  will  describe  the 
major  parts  of  NONISOCODE  and  discuss  special  features  in 
the  code  developed  for  this  study. 

Program  Overview 

NONISOCODE  is  a  two-dimensional  Navier-Stokes  simula¬ 
tion  designed  for  use  with  super-  or  hypersonic  flows.  It 
uses  MacCormack’s  explicit  method  to  compute  flow  quantities 
and  includes  a  pressure  damping  subroutine  to  allow  damping 
of  Gibbs  phenomena.  The  code  from  which  NONISOCODE  was 
adapted  was  originally  designed  to  analyze  three-dimensional 
flow  problems;  the  basic  structure  of  NONISOCODE  still 


reflects  this.  The  original  program  was  designed  to  compute 
flow  conditions  for  an  axial  'page*  and  then  sum  all  the 
'pages'  together  for  a  three-dimensional  flow  solution. 
This  structure  facilitated  vector i za t i on  of  the  code  for 
parallel  processor  computers;  however,  it  does  run  on  ’sca¬ 
lar'  machines  as  well.  NONISOCODE  only  computes  two-dimen¬ 
sional  flows,  but  still  retains  the  three-dimensional  struc¬ 
ture;  thus,  flow  variables  will  be  written  as  P(K,J,1).  The 
inclusion  of  the  unitary  third  dimension  reflects  the  actual 
two-dimensional  nature  of  the  problem.  A  brief  discussion 
of  each  subroutine  follows. 

Subroutine  MAIN 

This  portion  of  NONISOCODE  is  the  executive.  It  ac¬ 
cepts  input  values,  directs  most  subroutine  calls,  and  out¬ 
puts  results.  There  are  three  output  files.  The  first  is 
used  to  restart  the  computation  from  a  given  point  in  time. 
This  is  done  because  a  converged  solution  generally  requires 
about  40,000  iterations;  most  computer  operating  systems 
have  much  lower  time  limits  than  implied  by  that  number. 
The  second  output  file  is  tied  to  this  restart  capability. 
In  order  to  get  a  quick  look  at  the  state  of  the  computed 
flow  field,  this  output  file  holds  the  most  important  flow 
variables  at  each  point  in  the  field;  density,  velocity 
components,  pressure,  temperature,  and  turbulent  eddy 
viscosity.  In  addition,  there  is  a  list  of  surface  condi- 


tions  for  each  grid  point  along  the  surface.  The  final  file 
is  a  graphical  output  file  which  is  designed  to  interface 
with  several  common  graphics  packages.  This  file  contains 
velocity  components,  temperature,  pressure  ratios,  Mach 
number,  and  density  for  each  grid  point  in  the  field. 

Subroutine  EDDY 

This  subroutine  computes  the  eddy  viscosity  matrix  for 
the  flow  field  using  the  Ba Id w i n- Lorn ax  viscosity  model. 
However,  since  wind  tunnel  tests  have  shown  that  the  flow 
over  the  wedges  modeled  is  laminar,  this  subroutine  is 
superfluous  and  is  bypassed.  Subroutine  MAIN  loads  zero 
values  into  the  eddy  viscosity  matrix  and  does  not  call 
EDDY. 

Subroutine  PREAMB 

PREAMB  is  used  to  establish  the  initial  conditions  for 
the  solution  of  the  Navier-Stokes  equations.  PREAMB  estab¬ 
lishes  the  freestream  flow  field  over  the  entire 
computational  domain,  thus  giving  an  impulsive  start  to  the 
solution.  As  will  be  discussed  later,  the  initial  condi¬ 
tions  at  the  upstream  boundaries  permit  the  use  of  a  wedge 
angle.  There  is  a  provision  for  three  points  in  the 
freestream  to  allow  the  fourth-order  damping  routine  to 
operate  in  the  stagnation  region.  The  surface  initial  con¬ 
ditions  are  no-slip.  The  downstream  boundary  is  nonreflec- 
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tive,  allowing  waves  to  pass  out. 


None  of  the  other 
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boundaries  are;  this  can  be  a  cause  of  many  problems  if 
grid  is  improperly  designed. 


the 


Subroutine  TMSTEP 

TMSTEP  computes  the  CFL  condition  required  for  stabili¬ 
ty  of  the  solution.  The  CFL  condition  is  computed  for  each 
space  step  in  a  streamwise  row  of  grid  points.  The  minimum 
time  step  in  each  row  is  selected  and  compared  to  other 
minimum  time  steps  in  the  grid  from  the  surface  to  the  far- 
field  boundary.  The  minimum  CFL  condition  is  then  chosen  as 
the  CFL  condition  for  the  entire  flow  field.  It  is  then 
multiplied  by  an  input  CFL  factor  (always  less  than  one)  to 
guarantee  numerical  stability.  A  CFL  factor  is  required 
because  the  CFL  condition  only  guarantees  neutral  numerical 
stability  and  because  the  diffusion  terms  weren't  included. 
From  this  description,  it  can  be  seen  that  such  a  scheme  is 
inefficient  with  grids  that  are  irregular  or  are  highly 
stretched.  The  minimum  CFL  condition  will  be  based  on  the 
smallest  space  step  but  will  drive  the  integration  step  size 
for  the  largest  space  step.  Consequently,  the  solution  will 
slow  down  proportionally.  This  was  the  case  in  this  study; 
the  grid  used  is  semi-adaptive  and  exponentially  stretched 
in  the  normal  direction.  The  CFL  condition  will  be  based  on 
space  step  sizes  near  the  wall,  while  space  steps  at  the  far 
field  boundary  are  an  order  of  magnitude  larger. 
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Subroutine  BC 


This  subroutine  sets  the  boundary  conditions  for  the 
solution.  It  is  called  for  each  iteration.  The  surface 
boundary  condition  is  the  classical  no-slip  boundary  condi¬ 
tion;  this  prevents  the  investigation  of  possible  slip  flow 
in  the  stagnation  region  of  the  wedge.  All  inflow  boundary 
conditions  are  Dirichlet  boundary  conditions  where  the 
values  of  the  flow  properties  are  specified.  The  outflow 
boundary  (downstream)  uses  Neumann  boundary  conditions, 
where  the  gradients  of  the  flow  properties  are  zero.  Thus, 
all  inflow  boundaries  do  not  permit  waves  to  reflect  from 
them  (which  isn't  physically  correct  and  can  cause  numerical 
problems);  having  a  shock  wave  impinge  on  the  far-field 
boundary,  for  instance,  is  a  sure  way  to  bomb  the  solution. 
The  downstream  boundary  does  permit  waves  to  be  absorbed  and 
not  reflected;  a  shock  can  'exit’  the  downstream  boundary 
and  not  affect  the  solution.  The  use  of  Dirichlet  boundary 
conditions  for  inflow  boundaries  (upstream  and  far-field) 
specify  the  minimum  height  of  the  normal  grid;  it  must  be 
higher  at  the  trailing  edge  than  any  waves  that  might  impinge. 

Subroutine  TRANS 

TRANS  is  used  to  transform  the  physical  coordinate 
system  into  a  regular,  rectangular  grid  for  computation.  A 
generalized  transformat  ion  is  derived  using  the  chain  rule 
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(assuming  that  a  transformation  exists,  i.e.,  the  Jacobian 
of  the  transformation  is  nonzero  and  positive).  The 
transformation  scale  factors  (derivatives)  are  approximated 
using  second  order  central  differences.  The  output  from 
TRANS  is  a  transformed  grid  matrix  of  scale  factors  used  in 
all  subsequent  computations. 

Subroutines  PAGE ,  LETA(J),  LZETA(J),  SUM 

This  set  of  subroutines  represents  the  application  of 
MacCormack's  explicit  method  to  the  solution  of  the  two- 
dimensional  Navier-Stokes  equations.  PAGE  is  a  holdover 
from  the  original  three-dimensional  version  of  the  code;  it 
is  essentially  a  controller  for  the  following  subroutines. 
LETA(J)  computes  the  flow  derivatives  with  respect  to  the 
streamwise  transformed  coordinate  eta;  LZETA(J)  does  the 
same  with  respect  to  the  tranformed  normal  coordinate  zeta. 
SUM ( J )  sums  both  of  those  sets  of  derivatives  for  a  final 
solution.  SUM(J)  also  calls  the  final  subroutine,  when 
appropriate . 

Subroutine  DAMPING 

DAMPING  is  used  to  control  the  generation  of  Gibb's 
phenomena  (oscillations)  in  the  solution.  The  type  of  pres¬ 
sure  damping  used  is  an  artificial  viscosity.  Although  the 
Mach  number  modeled  is  quite  high,  the  pressure  damping 
required  turned  out  to  be  rather  low,  on  the  order  of  0.4. 


28 


SPECIAL  FEATURES 


Two  special  features  were  incorporated  into  NONISOCODE 
in  order  to  analyze  the  non i sothe rma  1  wall.  The  first  is 
the  incorporation  of  the  wedge  angle.  The  original  code 
from  which  NONISOCODE  was  adapted  could  only  consider 
axisyrametric,  zero  angle  flows.  In  order  to  perform  this 
study,  wedge  angles  had  to  be  included  in  subroutines  PREAMB 
and  BC.  In  those  subroutines,  the  velocity  components  U  and 
V  were  replaced  by  U  cos  S  and  -U  sin  &.  Once  the  wedge 
angle  is  included  in  PREAMB  and  BC,  the  flow  solution  will 
include  it  from  then  on. 

The  other  special  feature  included  in  NONISOCODE  was 
the  ability  to  enter  a  step  in  wall  temperature  at  any  point 
along  the  surface  of  the  wedge.  This  required  two  modifica¬ 
tions.  The  first  was  a  major  change  in  some  of  the  logic  in 
PREAMB  and  BC.  The  step  in  wall  temperature  was  included 
with  IF  statements  that  keyed  on  where  along  the  surface  the 
solution  was.  Before  a  certain  point,  Tw  would  be  Tw^; 
after  that  point,  Tw2*  The  second  modification  installed  a 
key,  KG,  which  corresponds  to  the  station  of  the  wall 
temperature  discontinuity,  in  the  input  files  and  the  code. 
Thus,  a  user  can  set  Twl,  Tw2,  and  KG  to  whatever  is  de¬ 
sired;  the  wall  temperatures  can  be  the  same  (as  was  done  in 
this  study  for  the  isothermal  wedge)  or,  for  the  same  ef¬ 
fect,  KG  can  specify  the  trailing  edge  of  the  surface.  Or, 
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as  was  done  for  the  nonisothermal  portion  of  the  study,  the 
discontinuity  can  be  placed  anywhere  on  the  surface.  In 
order  to  preserve  the  generality  of  the  code,  Twl  and  Tw2 
are  specified  in  separate  statements  that  determine  Tw, 
which  is  then  used. 

STEADY 

STEADY  is  the  heat  transfer  code  developed  specifically 
for  this  study.  It  is  a  simple  code,  but  it  does 
incorporate  some  interesting  points.  A  listing  of  STEADY 
can  be  found  in  Appendix  A. 

STEADY  is  designed  to  read  the  restart  files  of  NONISO¬ 
CODE  and  extract  the  field  temperature  from  them.  It  then 
computes  the  reference  convective  heat  transfer  coefficient 
to  be  used  to  nondimensionalize  the  convective  heat  transfer 
coefficient.  The  convective  heat  transfer  coefficient  is 
then  computed.  Again,  the  change  in  wall  temperature  is 
accounted  for  (if  the  wedge  is  nonisothermal;  if  not,  an 
input  value  bypasses  those  lines).  The  coefficient  of  vis¬ 
cosity  is  computed  based  on  the  wall  temperature.  Thermal 
conductivity  thus  includes  the  varying  wall  temperature. 

One  particularly  interesting  aspect  of  STEADY  is  the 
way  the  partial  derivative  of  temperature  with  respect  to 
normal  distance  is  computed.  Ty  must  be  computed  in  the 
transformed  plane  to  maintain  second  order  accuracy.  As 
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noted  before,  Ty  can  be  represented  as: 
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Due  to  the  physical  grid  used,  =  0.  Thus,  Ty  is  now 


represented  as: 
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This  derivative  can  be  computed  as  the  quotient  of  two 
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three-point,  one-sided  differences  with  respect  to  eta. 
Those  differences  are  second  order  accurate. 

There  are  two  output  files.  The  first  prints  a  listing 
of  X/L,  reference  convective  heat  transfer  coefficient, 
local  convective  heat  transfer  coefficient,  and  the 
nondimensional  ratio  H/Href.  The  second  output  file  is  a 
listing  of  X/L  and  H/Hreg  for  use  in  plotting. 

Convergence  is  determined  by  comparing  the  convective 
heat  transfer  coefficient  between  NONISOCODE  runs,  which  are 
usually  2000  iterations  apart.  If  the  percent  difference 
between  the  two  convective  heat  transfer  coefficients  is 
less  than  2%,  the  solution  is  considered  to  have  converged. 
The  difference  between  isothermal  and  nonisothermal  convec¬ 
tive  heat  transfer  coefficients  is  determined  in  essentially 
the  same  fashion.  The  equation  used  is: 


/  ^new  ^old  . 
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CHAPTER  V:  RESULTS 


OVERVIEW 

This  chapter  will  discuss  the  results  of  the  numerical 
simulation  of  the  isothermal  and  nonisothermal  wedges. 
There  will  be  a  brief  discussion  of  the  limitations  of  the 
results,  due  primarily  to  the  analytical  technique  used. 
The  heat  transfer  solution  for  both  the  isothermal  and  the 
nonisothermal  surfaces  will  be  discussed.  The  flow  field 
conditions  at  convergence  will  be  described.  The  presence  of 
shock  wave-boundary  layer  interaction  will  be  discussed 
briefly,  also.  Its  possible  effect  on  the  solution  will  be 
considered.  Finally  some  peculiarities  of  the  NONISOCODE 
solution  will  be  discussed 

FLOW  CONDITIONS 

Both  the  isothermal  and  nonisothermal  wedges  were 
modeled  under  the  same  conditions.  The  freestream  Mach 
number  was  14.24  and  the  freestream  stagnation  temperature 
was  2000  deg  R.  As  a  result,  the  freestream  static  tempera¬ 
ture  is  only  48.13  deg  R  and  the  freestream  static  pressure 
is  only  0.4981  lbs/sq  ft.  The  speed  of  sound  is  340  ft/sec. 
The  wedges  were  modeled  with  a  three  degree  deflection 
angle,  the  lowest  angle  discussed  in  Cappelano  (1:14). 
These  conditions  guarantee  that  boundary  layer  gradients 

will  be  quite  steep  and  the  stagnation  area  conditions 
severe.  The  surface  temperatures  were  as  follows.  For  the 


isothermal  wedge,  the  surface  was  considered  to  be  at  a 
uniform  temperature  of  540  deg  R.  The  nonisothermal  wedge 
had  a  leading  edge  temperature  of  540  deg  R  and  a  trailing 
edge  temperature  (after  the  discontinuity)  of  700  deg  R. 

LIMITATIONS  OF  THE  RESULTS 

The  chief  limitation  in  this  study  is  the  size  of  the 
streamwise  integration  step,  delta  x.  Considerable  time  was 
spent  on  determining  and  applying  the  correct  distribution 
of  normal  grid  points.  This  was  necessary  in  order  to 
sufficiently  resolve  the  temperature  inversion  layer  in  the 
region  of  interest  near  the  discontinuity.  In  the  course 
of  that  investigation,  it  was  found  that  NONISOCODE  did  not 
handle  a  highly  stretched  grid  well;  in  fact,  if  the 
exponential  stretch  factor,  Q,  exceeded  about  -0.15,  this 
particular  problem  wouldn't  run  at  all.  That  is  due  to  a 
combination  of  the  grid  and  the  particular  initial  condi¬ 
tions  (high  Mach  number,  steep  gradients,  etc.).  Due  to  the 
time  spent  on  the  normal  grid,  there  was  insufficient  time 
to  properly  investigate  the  streamwise  grid  distribution. 
Because  a  discontinuity  in  surface  temperature  is  being 
investigated,  it  would  be  desirable  to  have  a  high  concen¬ 
tration  of  grid  points  in  the  vicinity  of  the  discontinuity. 

In  fact,  at  one  point,  a  double  exponential  stretch  was 
considered,  with  the  first  stretch  in  the  normal  direction 


and  the  second  extending  both  up-  and  downstream  from  the 


discontinuity.  However,  the  limitations  of  NONISOCODE  and 
the  lack  of  time  available  forced  the  use  of  a  constant  grid 
distribution  in  the  streamwise  direction.  Thus,  the  solu¬ 
tion  in  the  immediate  vicinity  of  the  discontinuity  suffers 
from  truncation  error  and  subsequent  'smearing*.  As  will  be 


noted  later,  this  effect  can  be  seen  in  the  plots  of  the 
nonisothermal  convective  heat  transfer  coefficient. 

Another  limitation  to  the  results  to  be  discussed  is 
the  modeling  of  the  leading  edge.  When  the  normal  grid 
distribution  was  established,  it  was  determined  that  a  semi- 
adaptive  grid  would  be  used,  one  that  grew  in  relation  to 
the  growth  of  the  boundary  layer.  It  was  desired  to  main¬ 
tain  about  ten  points  in  the  boundary  layer  to  sufficiently 
resolve  the  temperature  inversion.  The  profile  that  the 
grid  was  based  on  was  at  X/L  =  0.2.  This  was  done  in  order 
to  escape  the  direct  effects  of  the  shock  wave-boundary 
layer  interaction  and  the  stagnation  region  and  to  avoid  a 
more  stringent  stability  criterion.  As  a  result,  from  X/L  = 
0  to  X/L  =  0.2,  the  normal  grid  distribution  is  constant  and 
will  not  necessarily  contain  enough  points  in  the  boundary 
layer  to  accurately  model  the  temperature  inversion.  When 
heat  transfer  is  computed,  the  partial  derivative  is 

computed  using  three  point,  one-sided  differences  in  ^ and 
.  This  requires  at  least  three  points  inside  the  tempe¬ 
rature  inversion  to  approximate  the  slope.  Those  points  are 


not  guaranteed  prior  to  X/L  =  0.2.  Thus,  no  results  are 
presented  for  the  leading  edge  region;  the  numbers  would  be 
inaccurate  and  misleading. 
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Another  point  to  consider  in  interpreting  the  results 
is  the  format  in  which  the  results  are  presented.  Strictly 
speaking,  this  isn't  a  limitation,  but  it  is  important  to 
correctly  interpreting  the  results.  Consider  Figures  4  and 
5;  they  are  representative  surface  and  contour  plots, 
respectively.  Note  that  in  Figure  4,  the  data  appears  to 
stop  along  a  parabolic  curve,  nowhere  near  the  apparent  far- 
field  boundary  shown  in  the  plot.  In  fact,  the  parabolic 
curve  _is  the  far-field  boundary  of  the  computational  domain, 
as  explained  earlier.  The  appearance  of  the  plot  is  an 
artifact  of  the  plotting  program,  DISPLAA.  DISPLAA  only 
recognizes  a  rectangular  grid,  so  that  data  presented  for 
the  parabolic  semi-adaptive  grid  shows  up  superimposed  on  a 
zero-value  rectangular  background. 

The  same  effect  is  visible  in  Figure  5,  the  representa¬ 
tive  contour  plot.  Here,  the  edge  of  the  computational 
domain  is  shown  by  a  very  close  clustering  of  contours, 
showing  the  drop  from  free-stream  values  to  an  artificial 
zero  value.  in  both  cases,  the  zero-value  data  had  to  be 
inserted  to  allow  DISPLAA  to  plot  the  semi-adaptive  grid 
results.  As  a  beneficial  side  effect,  the  presence  of  the 
zero-value  data  serves  to  constantly  highlight  the  shape  of 
the  physical  computational  domain. 


HEAT  TRANSFER  SOLUTIONS 


The  first  results  to  be  discussed  are  the  heat  transfer 
solutions.  Figure  3  is  a  plot  of  the  nond i mens i ona 1 i zed 
convective  heat  transfer  coefficients  for  the  isothermal  and 
the  non  isothermal  surfaces  and  experimental  data  from  Hodge 
et  al  (2:8,  10).  Nond imens iona 1 i zat ion  was  accomplished  by 
ratioing  the  dimensional  convective  heat  transfer  coeffi¬ 
cients  to  a  reference  convective  heat  transfer  coefficient 
described  earlier.  Hodge's  data  was  already 
nondimensionalized  in  the  same  fashion;  the  only  correction 
that  needed  to  be  made  was  to  shift  his  data  points 
downstream  to  account  for  the  presence  of  three  grid  points 
in  the  freestream  ahead  of  the  leading  edge  in  this  study. 

Consider  first  the  isothermal  convective  heat  transfer 
coefficient,  plotted  as  a  solid  line  in  Figure  3.  It  is 
fairly  constant  about  the  value  h/hre£  =  1.741  for  the 
region  from  X/L  =  0.2  to  X/L  =  0.9.  There  is  only  a  3.4% 
variation  from  this  value  over  this  range.  Note  that  the 
ratio  does  tend  to  trail  off  somewhat  towards  the  trailing 
edge  of  the  wedge.  The  theoretical  value  for  an  isothermal 
h/href  is  constant;  the  variation  shown  is  probably  due  to 
truncation  error  in  the  computations. 

The  nonisothermal  h/hre£  is  shown  by  the  dotted  line  in 
Figure  3.  Up  to  about  X/L  =  0.42,  the  nonisothermal  and  the 
isothermal  values  are  identical,  which  is  to  be  expected. 
At  X/L  =  0.42,  the  nonisothermal  curve  begins  to  oscillate 
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somewhat  in  the  immediate  vicinity  of  the  temperature 
discontinuity  on  the  wedge.  Immediately  thereafter,  h/href 
plummets  to  0.64  and  then  recovers  slowly  to  values  about 
17%  below  the  isothermal  values.  This  difference  between 
the  isothermal  and  the  non i sot  he rma 1  values  after  the 
d  iscont i nu i ty  is  expected  due  to  the  higher  temperature  of 
the  trailing  edge's  surface,  160  deg  R  hotter  than  for  the 
isothermal  wall. 

The  slow  recovery  of  the  nonisothermal  h/hre^  is  due  in 
part  to  truncation  error  'smearing'  the  convective  heat 
transfer  coefficient.  However,  some  of  the  lag  in  recovery 
is  due  to  the  axial  diffusion  of  energy.  Within  the  scope 
of  this  study,  it  is  impossible  to  separate  the  contribu¬ 
tions  of  truncation  error  and  axial  diffusion  to  the  h/hre£ 
recovery  lag;  it  is,  however,  clear  that  axial  diffusion  of 
energy  is  occuring,  that  the  nonisothermal  wall  effect  is 
present.  The  lack  of  recovery  indicates  that  some  of  the 
energy  is  being  convected  downstream  and  is  consequently 
reducing  the  local  convective  heat  trarsfer. 

Hodge's  experimental  data  are  also  presented  in  Fig.  3. 
The  isothermal  data  is  for  a  thin  stainless  steel  plate 
(2:2,  5);  the  nonisothermal  data  is  for  an  FRSI  test  sample 
(2:2,  6).  The  isothermal  data,  plotted  as  +'s,  shows  a 
close  correlation  with  the  isothermal  Navier-Stokes  predic¬ 
tion,  with  only  a  2.3%  average  variation  from  the  computed 
data  over  the  trailing  edge  portion  of  the  plate.  The 


agreement  is  best  for  the  first  five  points;  as  the  trailing 
edge  is  approached,  the  experimental  data  diverges  slightly. 
The  reason  for  this  divergence  is  unknown.  There  are  only 
three  data  points  available  for  the  non  isothermal  FRSI  test 
article.  Fig.  3  shows  these  points,  plotted  as  0's, 
reflecting  the  same  trend  of  recovery  after  the  discontinui¬ 
ty  as  the  Nav ier-Stokes  prediction.  The  magnitudes  of  the 
computed  and  experimental  convective  heat  transfer 
coefficients  match  fairly  well  for  the  first  and  third 
experimental  points;  the  middle  experimental  point  is  well 
off  the  computed  curve.  The  reason  for  this  disagreement  is 
unknown.  Note,  however,  that  the  behavior  of  the  two 
computed  curves  and  the  isothermal  thin-skin  points  is 
similar;  they  all  seem  to  have  a  "hump"  between  X/L  =  0.6  and 
X/L  =  0.75. 

Finally,  note  that  the  behavior  of  the  computed  and 
experimental  curves  is  quite  different.  The  computed 
isothermal  and  non  isothermal  curves  tend  to  follow  on 
another  after  the  discontinuity,  allowing  for  the  17% 
difference.  The  slope  and  shape  of  the  two  curves  are  quite 
similar.  The  isothermal  and  nonisotherma  1  experimental 
curves  do  not  tend  to  mirror  one  another.  As  noted  above, 
the  nonisothermal  data  tends  to  approach  the  isothermal  data 
as  the  trailing  edge  of  the  wedge  is  approached.  indeed, 
the  last  values  show  only  a  7.6%  difference  between  them.  It 
has  been  suggested,  in  another  context,  that  what  is 


happening  in  the  experimental  data  near  the  trailing  edge  is 
a  three-dimensional  effect.  While  there  is  very  little 
concrete  evidence  to  go  on  in  this  study,  it  is  safe  to  say 
that  the  experimental  and  computed  non i sothermal  values  near 
the  trailing  edge  are  in  fairly  close  agreement.  This  would 
seem  to  mitigate  against  any  three-dimensional  effects. 

FLOW  FIELD  SOLUTIONS 

Results  of  NONISOCODE  simulations  are  plotted  in  Fi¬ 
gures  4  through  27.  The  results  for  the  isothermal  wall  are 
shown  in  Figures  4  through  15;  the  nonisothermal  wall  in 
Figures  16  through  27.  The  surface  plots  for  pressure  and 
density  clearly  show  the  leading  edge  shock  wave.  The  shock 
angle  is  in  the  vicinity  of  6  degrees,  agreeing  with  theory. 
The  pressures  behind  the  shock  are  about  50%  higher  than 
predicted  by  Rankine-Hugoniot  shock  jump  equations;  this  is 
due  to  the  shock  wave-boundary  layer  interaction.  Behind 
the  shock,  outside  of  the  boundary  layer,  the  freestream 
static  temperature  is  close  to  the  64.95  deg  R  predicted  by 
Rankine-Hugoniot  theory.  As  these  examples  show,  NONISOCODE 
has  done  a  good  job  of  modeling  the  flow  field.  One 
interesting  point  is  the  size  of  the  temperature  inversion. 
The  surface  plots  of  temperature.  Figures  10  and  22,  show 
the  inversion  clearly  defined  by  the  ridge  near  the  bottom 
edge  of  the  plot.  It  can  be  seen  that  the  inversion  is  only 
about  0.0096  ft  tall  (about  0.1  in).  This  is  a  very  small 


region  for  such  a  large  temperature  rise. 

Another  interesting  point  is  the  surface  pressure  dis¬ 
tribution.  Figure  28  shows  the  isothermal  and  nonisothermal 
surface  pressure  distributions,  as  well  as  a  theoretical 
Rankine-Hugoniot  line.  The  analytical  curves  are  identical 
prior  to  X/L=0.6,  but  after  the  flow  crosses  the  thermal 
discontinuity,  the  nonisothermal  surface  pressure 
experiences  a  pressure  spike  of  about  0.2  lbs/sq  ft.  The 
spike  occurs  about  3/4  inch  (0.0632  ft)  behind  the 
discontinuity.  After  this  spike,  the  nonisothermal  surface 
pressure  decays  in  a  similar  fashion  to  the  isothermal 
surface  pressure,  but  at  a  value  about  4%  above  it. 

The  pressure  spike  can  be  explained  in  terms  of  the 
thermal  boundary  layer.  At  the  discontinuity,  the  surface 
temperature  jumps  from  540  deg  R  to  700  deg  R;  this  will 
cause  a  corresponding  increase  in  the  thickness  of  the 
thermal  boundary  layer.  When  this  "bloom"  occurs,  compres¬ 
sion  waves  are  generated  in  the  flow  field.  Assuming  that 
the  Mach  cone  for  these  compression  waves  is  centered  at  the 
edge  of  the  boundary  layer  at  the  discontinuity,  the  Mach 
angle  (defined  as  the  inverse  sine  of  the  reciprocal  of  the 
edge  Mach  number)  is  a  little  greater  than  4  degrees.  The 
compression  wave  will  impinge  the  surface  of  the  wedge  about 
0.6  inches  behind  the  discontinuity.  Referring  to  Figure 
28,  it  can  be  seen  that  this  is  where  the  pressure  spike 
occurs.  The  actual  calculated  point  of  impingment  is  at  X/L 


=  0.69564;  the  computed  pressure  spike  occurs  at  X/L  - 
0.6681.  This  4%  discrepancy  can  be  explained  by  the  fact 
that,  as  the  compression  wave  goes  deeper  into  the  boundary 
layer,  flow  Mach  number  decrease  rapidly  and  the  compression 
wave  is  "bent  back",  causing  the  wave  to  impinge  slightly 
ahead  of  the  predicted  point. 

The  Rank i ne-Hugon  iot  line  at  the  bottom  of  the  plot 
shows  the  theoretical  downstream  pressure;  if  traditional 
boundary  layer  assumptions  had  been  applied,  this  v  value 
would  also  represent  the  surface  pressure  behind  the  shock. 
It  serves  to  demonstrate  the  effect  of  the  shock  wave¬ 
boundary  layer  interaction;  the  predicted  pressure  behind 
the  shock  is  significantly  lower  than  the  computed  (and 
actual)  pressure  due  to  the  interaction  of  the  shock  wave 
and  boundary  layer. 

In  general,  the  flow  field  solutions  for  the  isothermal 
and  nonisothermal  wedges  are  quite  comparable,  with  the 
obvious  difference  of  the  temperature  discontinuity.  Both 
solutions  are  accurate  models  of  the  real  flow  field.  With 
the  agreement  between  model  and  reality  demonstrated,  the 
heat  transfer  results  can  be  considered  as  accurately 
modeling  the  physical  situation. 

SHOCK  WAVE- BOUNDARY  LAYER  INTERACTIONS 

An  interaction  of  the  bow  or  leading  edge  shock  wave 
and  the  boundary  layer  is  a  phenomenon  that  will  occur  in 


most  hypersonic  aerodynamics  problems.  Its  presence  in  this 
problem  was  not  unexpected.  The  presence  of  the  interaction 
affects  the  flow  field  solution  and  the  heat  transfer.  The 
discussion  that  follows  is  based  on  the  work  of  Dorrance 
(11:144-148).  It  should  be  noted  that  no  effort  was  made  to 
analyze  in  depth  the  shock  wave-boundary  interaction; 
rather,  the  presence  of  the  interaction  was  observed,  its 
strength  calculated,  and  a  qualitative  overview  of  its  ef¬ 
fects  was  made. 

When  an  object  is  moving  at  hypersonic  speeds,  the 
displacement  thickness  is  much  larger  than  for  the  same 
object  at  slower  speeds.  This  is  due  to  the  higher  tempera¬ 
tures  found  in  a  hypersonic  boundary  layer.  The  larger  size 
of  the  displacement  thickness  means  that  the  effective  shape 
of  the  body  is  changed  to  include  the  displacement 
thickness.  At  the  leading  edge,  the  flow  is  turned  sharply 
away  from  the  surface;  this  turning  gives  rise  to  a  series 
of  compression  waves  that  coalesce  into  a  shock.  The  very 
presence  of  the  shock  affects  the  growth  of  the  boundary 


layer.  This  interaction  is  the  shock  wave  boundary  layer 
interaction . 

The  strength  of  this  interaction  can  be  estimated  by 
computing  an  interaction  factor  %  ,  defined  as: 


where 


Y  =  (C1/2  M3)/Re1//2 
c  =  /7^/(  Pe/Ue) 


(5.2) 


If  this  interaction  factor  is  greater  than  1,  then  the 
effects  of  shock  wave-boundary  layer  interaction  must  be 
accounted  for.  In  the  current  problem,  at  the  leading 
edge  has  a  value  of  about  20.2;  it  drops  off  to  a  value  of 
2.65  at  the  trailing  edge  of  the  wedge.  ^  is  greater  than  1 
and  the  interaction  must  be  considered.  Referring  again  to 
Fig.  28,  the  presence  and  magnitude  of  the  interaction  can 
be  seen  by  comparing  the  analytical  pressure  curves  with  the 
theoretical  line.  Near  the  leading  edge,  the  computed  pres¬ 
sure  is  about  50%  greater  than  the  theoretical  value.  Fur¬ 
ther  downstream,  this  difference  decreases  to  about  20%. 

NUMERICAL  OSCILLATIONS 

When  examining  the  surface  plots  in  Figures  4  through 
27,  it  is  very  easy  to  see  oscillations  in  the  plots  for 
streamwise  velocity,  temperature,  and,  to  a  lesser  extent, 
Mach  number  and  density.  These  oscillations,  seen  most 
easily  in  the  surface  plots  for  streamwise  velocity,  have 
one  probable  origin. 

The  most  probable  explanation  is  that  the  oscillations 
are  artifacts  of  the  plotting  process.  This  determination 
comes  from  the  nature  of  DISPLAA,  the  program  used  to  gene¬ 
rate  the  plots.  The  version  of  DISPLAA  used  requires  a 
rectangular  grid  to  plots  its  data.  If  the  data  are 
presented  in  a  non-rectangular  format,  DISPLAA  uses  a  nine- 


point  linear  average  interpolating  scheme  to  interpolate  the 
data  to  a  rectangular  grid.  This  problem  did  not  use  a 
rectangular  grid,  and  the  results  are  presented  in  the  semi- 
adaptive  grid  format.  Thus,  DISPLAA  has  interpolated  the 
data  to  fit  its  own  grid.  Now  the  oscillations  are  most 
apparent  in  the  streamwise  velocity  surface  plots  and  the 
temperature  surface  plots.  The  gradients  here  from  the 
surface  to  the  freestream  are,  in  absolute  terms,  the  most 
severe  encountered  in  the  entire  solution.  The  next  most 
severe  gradients  are  to  be  found  in  the  temperature  solu¬ 
tion,  partly  because  temperature  is  dependent  on  flow  velo¬ 
city.  Realizing  that  the  data  is  interpolated  and  realizing 
that  the  most  severe  gradients  are  present  in  the  plots 
where  the  oscillations  are  most  apparent  lead  directly  to 
the  argument  that  those  oscillations  are  artifacts  of  DIS¬ 
PLAA  and  are  enhanced  by  steep  gradients. 

A  point  to  buttress  this  argument  is  that  the  oscilla¬ 
tions  do  not  show  up  in  the  data.  An  investigation  of  the 
raw  data  plotted  shows  no  oscillations  of  the  kind  seen  in 
the  surface  plots.  Also,  the  oscillations  do  not  show  up  in 
the  contour  plots.  Thus,  it  is  improbable  that  the  oscilla¬ 
tions  have  any  other  cause. 
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CHAPTER  VI:  CONCLUSIONS  AND  RECOMMENDATIONS 


OVERVIEW 

This  chapter  will  discuss  conclusions  and  recommenda¬ 
tions  derived  derived  from  the  study.  Conclusions  regarding 
the  numerical  methods  used  and  concerning  the  nonisothermal 
wall  effect  will  be  discussed.  Recommendations  for  further 
study  and  possible  methods  will  be  presented. 

CONCLUSIONS 


NONISOTHERMAL  WALL  EFFECT 

This  study  has  demonstrated  that  the  nonisothermal  wall 
effect  can  be  modeled  using  a  two-dimensional  Navier-Stokes 
simulation.  The  nonisothermal  wall  effect  manifests  itself 
as  an  'incomplete'  recovery  in  the  nonisothermal  convective 
heat  transfer  coefficient  after  the  temperature  discontinui¬ 
ty  on  the  wedge  surface.  The  actual  magnitude  of  the  noni¬ 
sothermal  wall  effect  was  not  determined  due  to  numerical 
'smearing'  caused  by  truncation  error;  however,  there  is  no 
doubt  that  the  nonisothermal  wall  effect  was  modeled  by 
NONISOCODE.  Additionally,  it  has  been  shown  that  the 
numerical  results  agree  closely  with  experimental  results; 
differences  can  be  explained  by  the  lack  of  axial  resolution 
in  the  region  of  the  thermal  discontinuity. 

A  subsidiary  conclusion  that  can  be  drawn  is  that 


NONISOCODE  adequately  models  hypersonic  flow  over  wedges. 
This  is,  of  course,  a  necessary  condition  for  the  proper 
simulation  of  the  non  isothermal  wall  effect.  The  tempera¬ 
ture  inversion  layer  and  the  axial  diffusion  of  heat  energy 
at  the  thermal  discontinuity  are  correctly  modeled  by 
NONISOCODE  . 

SHOCK  WAVE-BOUNDARY  LAYER  INTERACTION 

It  has  been  shown  that  a  shock  wave-boundary  layer 
interaction  did  occur  and  did  have  observable  effects  on  the 
computed  solution.  The  20%  to  50%  increase  in  surface 
pressure  over  theoretical  predictions  is  in  general  agree¬ 
ment  with  Nagamatsu  and  Sheer. 

PRESSURE  SPIKE  DUE  TO  NONISOTHERMAL  WALL 

The  observed  pressure  spike  is  a  real  phenomenon  that 
is  probably  a  result  of  the  nonisothermal  wall  effect.  Its 
probable  cause  is  the  impingement  of  a  compression  wave  on 
the  surface  of  the  wedge.  The  compression  wave,  in  turn,  is 
probably  generated  by  the  sudden,  almost  discontinuous 
growth  of  the  momentum  and  thermal  boundary  layers  at  the 
discontinuity.  This  pressure  spike  may  be  a  distinguihing 
characteristic  of  the  nonisothermal  wall  effect,  and  might 
be  used  to  signal  its  existence. 

NUMERICAL  METHODS 

The  numerical  methods  used  in  this  study  were  adequate 


to  properly  model  the  non  isothermal  wall  effect.  MacCor- 
mack's  explicit  method  was  appropriate  to  the  problem.  The 
grid  chosen  was  appropriate  to  the  problem,  with  two 
specific  exceptions  (to  be  addressed  below).  In  particular, 
the  use  of  a  semi-adaptive  grid  that  grows  (in  the  normal 
direction)  with  the  boundary  layer  was  essential  to  the 
success  of  the  study.  This  grid  permitted  the  investigation 
of  the  temperature  inversion  in  the  boundary  layer  while 
allowing  shock  and  Mach  waves  to  escape  the  downstream 
boundary  without  excessive  stretching. 

The  method  used  to  compute  the  convective  heat  transfer 
coefficient  was  entirely  appropriate;  the  use  of  three- 
point,  one-sided  finite  differences  in  the  transformed  plane 
for  <JT/3y  retained  second  order  accuracy.  The  method  used  to 
determine  convergence  is  simple  and  straight-forward,  but  it 
adequately  demonstrates  convergence  of  temperature,  the 
figure  of  merit  for  this  study. 

There  were,  however,  some  significant  problems  with  the 
numerical  methods  used.  The  most  significant  problem  was 
the  choice  of  grid.  The  use  of  a  regularly-distributed  grid 
in  the  streamwise  direction  simplified  computations  but  it 
also  induced  numerical  error  in  the  stagnation  region  and  in 
the  vicinity  of  the  temperature  discontinuity.  The  solution 
in  the  stagnation  region  suffered  because  the  grid  did  not 
allow  resolution  of  the  stagnation  region.  A  similar  prob¬ 
lem  occurred  at  the  discontinuity.  The  problem  is  more 
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serious  here  because  the  lack  of  resolution  coupled  with  the 
existence  of  the  discontinuity  induces  a  significant  trunca¬ 
tion  error.  Thus,  the  less-than-total  convective  heat 
transfer  coefficient  recovery  is,  in  part,  an  artifact  of 
the  truncation  error.  While  the  recovery  does  demonstrate 
the  non  isothermal  wall  effect,  it  is  impossible  to  determine 
the  magnitude  of  the  effect  in  this  study. 

Another  problem  that  is  of  lesser  significance  is  the 
exclusive  use  of  the  three-point,  one-sided  finite 
difference  to  model  the  partial  derivative  of  temperature 
with  respect  to  normal  distance.  While  this  finite 
difference  worked  quite  well  for  X/L  >  0.4,  it  is  unreliable 
for  the  area  near  the  leading  edge.  This  is  due  to  the 
decreasing  number  of  points  in  the  temperature  inversion. 
It  is  possible  that  using  a  two-point,  one-sided  finite 
difference  near  the  stagnation  region  would  yield  valid 
results,  although  at  the  cost  of  numerical  accuracy.  Exclu¬ 
sively  using  the  three-point,  one-sided  finite  difference 
denied  investigation  of  the  convective  heat  transfer  coeffi¬ 
cient  near  the  leading  edge. 

RECOMMENDATIONS 

There  are  several  recommendations  that  result  from  this 
study.  The  first  is  that  the  non  isothermal  wall  effect 
needs  more  study.  This  effort  has  shown  that  it  is  possible 
to  accurately  model  the  non i sothe rma  1  wall  effect;  it  has 


also  demonstrated  that  there  are  several  shortcomings  in  the 
methods  used.  More  study  needs  to  be  done  in  order  to 
correctly  model  all  aspects  of  the  nonisothermal  wall  effect 
and  to  determine  its  magnitude. 

Closely  tied  to  the  first  recommendation  is  the  second, 
which  is  the  use  of  a  more  appropriate  grid  for  modeling  the 
nonisothermal  wall  effect.  Although  this  study  showed  that 
the  inviscid  portion  of  the  hypersonic  flow  can  be  modeled 
with  a  fairly  coarse  grid,  the  stagnation  region  and  the 
temperature  discontinuity  require  a  much  finer,  specialized 
grid.  In  particular,  the  grid  should  employ  dual 
stretching:  the  normal  distribution  must  be  stretched  in 
order  to  include  the  temperature  inversion,  while  the 
streamwise  distribution  must  stretch  at  the  stagnation  point 
and  at  the  discontinuity.  The  stretching  at  the  disconti¬ 
nuity  should  be  both  up-  and  downstream  to  capture  the  true 
character  of  the  nonisothermal  wall  effect. 

The  use  of  such  a  grid  will  raise  problems  with  NONISO- 
CODE.  During  the  study,  several  variants  of  stretched  grids 
were  tried;  very  few  worked.  The  reason  for  all  the 
failures  is  unknown;  however,  it  became  apparent  that  NONI¬ 
SOCODE  in  its  present  form  cannot  use  highly  stretched  grids 
(with  a  stretching  exponent  much  greater  than  0  =  -0.15). 
Whether  this  is  due  to  the  transformation  from  the  physical 
to  the  transform  plane  or  to  the  implementation  of  MacCor- 
mack's  method  or  to  some  other  source  is  unknown.  That 
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aspect  requires  more  investigation. 

A  further  recommendation  is  to  further  investigate 
numerically  the  shock  wave-boundary  layer  interaction  and 
how  it  affects  the  nonisothermal  wall  effect.  The  interac¬ 
tion  has  a  tendency  to  increase  the  surface  pressure 
downstream  of  the  shock  and  can  influence  the  heat  transfer 
characteristics  of  the  flow.  The  grid  used  did  not  allow 
sufficient  resolution  of  the  stagnation  region  or  the  region 
about  the  thermal  discontinuity.  A  different  grid  will  need 
to  be  used  to  investigate  this  problem. 

A  final  recommendation  is  to  further  study  the  pressure 
spike  that  is  seen  to  occur  downstream  of  the  discontinuity. 
It  is  possible  that  this  spike  is  a  distinguishing 
characteristic  of  the  nonisothermal  wall  effect.  Further 
investigation  could  show  whether  or  not  the  pressure  spike 
can  signal  the  presence  of  a  nonisothermal  wall  effect. 
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APPENDIX  A:  PROGRAM  LISTINGS  FOR  NONISOCODE  AND  STEADY 


PROGRAM  MAIN  ( INPUT , OUTPUT , TAPE5= INPUT ,TAPE6  =  OUTPUT , TAPE  1 ,TAPE2 
1  TAPE3 , TAPE4 ) 

THIS  IS  NONISOCODE,  A  2-D  NAVI ER-STOKES  SIMULATION. 
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CONTROL  THE  DATA  FLOW  AND  DIFFERENCE  OPERATORS 
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COMMON  /DEP/  RHO (62,30,1) , RHOU (62,30,1) , RHOV (62,3 
RHOW( 1,1,1) , RHOE( 62,30,1) 

DIMENSION  UET (62) ,UZT(62) ,U(62),V(62),C(62) 


DTC ( 1 )  =  1 .0 
GAMM2=GAMMA*GAMM1 
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COMPARING  DTMIN  BETWEEN  ADJACENT  PLANES 
DTCFL=AMINl ( DTCFL, DTMIN) 
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DO  1  K=2 , KLM 
IF  (K.LE.KG)  TW=TW 
IF  (K.GT.KG)  TW=TVi 
RH= 1 . O/RHOP { K , 2 , 1 ) 
RHOUP ( K , 1 ,1 )=0.0 


RHOVP { K , 1 ,1 )=0.0 

RHOP (  K  ,  1  , 1 )  = RHOP (  K  ,  2 , 1 )  *  ( RHOEP (  K  ,  2 , 1  )  *  RH-  0 . 5  * 

(  ( RHOUP ( K  f  2  r 1 ) *RH ) *  *  2+ ( RHOVP (K,2,l)*RH)**2))/( CV*  TW ) 
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DO  3  J=1,JL 

RHOP  ( 1 , J, 1 ) =RHOINF 

RHOUP { 1 , J , 1 ) =RUINF*COS( DELTA) 

RHOVP { 1 , J , 1 ) =-RUINF*SIN ( DELTA) 

RHOEP ( 1 , J , 1 ) =REINF 
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DO  14  K=2,KF 

RHO{ K , 1 , 1 )=RHOINF 

RHOU (  K  ,  1  , 1 )  =RUINF*COS( DELTA) 

RHOV{ K ,1 ,1 ) =-RUINF*SIN{ DoLTA) 

RHOE ( K  ,  1  , 1  ) =REINF 
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DO  103  JV=2,JLM 
JP=JV+1 
JM=JV- 1 

DO  104  KV=1,KL 

XET(KV,JV)=(X( KV , JP , 1 ) -X ( KV,  JM , 1 ) ) *RDET 


YET ( KV , JV ) = ( Y ( KV /JP,1)-Y( KV ,  JM , 1 ) ) *RDET 
104  CONTINUE 
103  CONTINUE 
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DO  301  KV=1,KL 

DJ(KV, JV)=XET(KV,JV)*YZT(KV,JV)-XZT(KV, JV)*YET(KV, JV) 
301  CONTINUE 
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COMMON  /TV/  ETX( 62,30,1) , ETY (62,30,1) , ZTX( 62,30,1) , ZTY (62,30,1) 
COMMON  /DEP/  RHO{ 62,30,1) , RHOU (62,30,1) , RHOV (62,30,1) , 


RHOW( 1,1,1) , RHOE (62,30,1) 

COMMON  /DEPP/  RHOP(62, 30, 1) ,RHOUP(62, 30,1) ,RHOVP(62, 30,1) 
RHOWP ( 1,1,1) ,RHOEP( 62,30,1) 


CALL  LZETA(J) 
CALL  SUM ( J ) 
CONTINUE 
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Wind  tunnel  tests  of  Space  Shuttle  Orbiter  insulating  articles  have  demonstrated 
presence  of  a  nonisothermal  wall  effect,  which  is  a  lag  in  heat  transfer  recovery  aftc.'*’ 
the  flow  passes  over  a  surface  temperature  discontinuity  resulting  in  a  downstream 
transport  of  energy.  Theoretical  analyses  and  numerical  simulations  of  hypersonic  flow 
over  discontinuous  nonisothermal  surfaces  using  boundary  layer  theory  have  also  indicated 
the  presence  of  this  effect. 

This  thesis  studies  the  nonisothermal  wall  effect  by  modeling  the  hypersonic  flow 
over  an  inclined  wedge  with  a  discontinuous  nonisothermal  surface.  The  flow  is  modeled 
using  the  two-dimensional  Navier-Stokes  equations.  MacCormack's  method  is  used  to  solve  ! 
the  Navier-Stokes  equations.  The  program  used  to  implement  these  methodologies  is 
discussed  and  a  listing  given.  A  semi -adapti ve  grid  is  used  to  represent  the  physical 
conditions  of  the  problem.  Heat  transfer  is  presented  as  a  nondimensional  ratio  of  the 
local  convective  heat  transfer  coefficient  to  a  reference  heat  transfer  coefficient. 

The  results  of  this  study  show  that  the  nonisothermal  wall  effect  can  be  successfully 
modeled  using  the  two-dimensional  Navier-Stokes  equations  and  MacCormack's  explicit 
method.  The  lag  in  the  recovery  of  the  convective  heat  transfer  coefficient  is  found  to 
match  lags  seen  in  other  analyses  of  the  problem.  Due  to  the  high  Mach  number  modeled,  a 
shock  wave-boundary  layer  interaction  is  found  to  have  an  effect  on  heat  transfer.  A 
significant  surface  pressure  spike  is  found  to  occur  downstream  of  the  discontinuity, 
which  may  be  a  result  of  the  increase  in  size  of  the  momentum  and  thermal  boundary  layers 
at  the  discontinuity. 

The  study  concludes  that  the  nonisothermal  wall  effect  can  be  adequately  modeled  by 
the  two-dimensional  Navier-Stokes  equations;  that  the  shock  wave-boundary  layer 
interaction  does  have  an  effect  on  the  heat  transfer;  and  that  the  occurence  of  a  spike  in 
surface  pressure  may  be  a  unique  result  of  the  nonisothermal  wall  effect.  Significant 
recomnendations  Include  the  need  for  further  study  of  the  nonisothermal  wall  effect,  the 

need  to  use  more  optimal  grids  and  solution  methods,  the  need  to  more  thorough _ 

Investigate  the  shock  wave-boundary  layer  effect,  and  the  need  for  further  study  of  tfie 
surface  pressure  response  to  the  nonisothermal  wall  effect. 
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